k=3;
%入射粒子初始位置
y0=[-10,1,10,0;
    -10,1,-10,0;
    -10,1,2,0;
    -10,1,-2,0;
    -10,2,2,0;
    -10,2,-2,0;
    -10,1,4,0;
    -10,1,-4,0;
    -10,1,15,0;
    -10,1,-15,0;
    -10,1,7,0;
    -10,2,-7,0;
    -10,2,7,0;
    -10,1,-7,0;
    -10,1,0,0.4;
    -10,1,0,-0.4;
    ];
%粒子运动方程
F=@(t,y)[y(2);
         k*y(1)/sqrt(y(1).*y(1)+y(3).*y(3)).^3;
         y(4);
         k*y(3)/sqrt(y(1).*y(1)+y(3).*y(3)).^3];
%画靶粒子
line(0,0,'marker','.','markersize',50,'color','k');
text(2,0,'重荷');
axis([-10 15 -20 20])
hold on
%粒子入射
for i=1:length(y0)
    [t,y]=ode45(F,[0:0.1:30],y0(i,:));
    plot(y(:,1),y(:,3));
end